###########################################
# Figure 2: Dynamic Average Causal Responses and least wiggly path of confound

###########################################

##
# The least winggly path of confound is computationally heavy. Thus, one can decide to produce the figure without it. 
with_LWP <- TRUE

##


# Define Change.Ded x year for parallel trend graphs
df$decrease_itm_2014 <- ifelse(df$year == 2014,df$itm_prop_change,0)
df$decrease_itm_2015 <- ifelse(df$year == 2015,df$itm_prop_change,0)
df$decrease_itm_2016 <- ifelse(df$year == 2016,df$itm_prop_change,0)
df$decrease_itm_2017 <- ifelse(df$year == 2017,df$itm_prop_change,0)
df$decrease_itm_2018 <- ifelse(df$year == 2018,df$itm_prop_change,0)
df$decrease_itm_2019 <- ifelse(df$year == 2019,df$itm_prop_change,0)
df$decrease_itm_2020 <- ifelse(df$year == 2020,df$itm_prop_change,0)
df$decrease_itm_2021 <- ifelse(df$year == 2021,df$itm_prop_change,0)
df$decrease_itm_2022 <- ifelse(df$year == 2022,df$itm_prop_change,0)


# Define main variables
dep <- "perc_yes ~"
Xs_elFE <-  paste("turnout","Bond","failed_recently",sep = "+")
FEs <- "|LEAID + election_FE|0"
cluster <- "|LEAID "
##

# Restrict to school district with more than 1 referendums (cannot identify with school district fixed effects)
dfreg <- df[df$one_election ==0,]
# Restrict to close elections as defined in Cellini et al. (2010)
dfreg25 <- dfreg[dfreg$tight_election ==1,]

## Event study 
main <- "  decrease_itm_2014+decrease_itm_2015 +decrease_itm_2016 + decrease_itm_2017 +
          decrease_itm_2019 +decrease_itm_2020 +decrease_itm_2021 +decrease_itm_2022 +"
R_pt <- felm(data = dfreg25,weights = dfreg25$population_votingage,formula = as.formula(paste(dep,main,Xs_elFE,FEs,cluster,sep = "")))

# Extract point estimates and confidence intervals
graphs <- matrix(nrow = 9,ncol = 4)
graphs[,2] <- c(R_pt$coefficients[1:4],0,R_pt$coefficients[5:8] )
graphs[,c(1)] <- c(confint(R_pt,level = 0.95)[1:4,1],0,confint(R_pt,level = 0.95)[5:8,1] )
graphs[,c(3)] <- c(confint(R_pt,level = 0.95)[1:4,2],0,confint(R_pt,level = 0.95)[5:8,2] )

# Organize dataframe for plotting in ggplot2
graphs <- as.data.frame(graphs)
graphs[,4] <- seq(2014,2022,1)
graphs$date <- as.Date(paste(graphs[,4],"01","01",sep = "-"), format = "%Y-%m-%d")

colnames(graphs) <- c("lb","coef","up","year","date")

### 
# Compute p-values of pre and post-trend test 
###

# Pretrend test
main_test <- " decrease_itm_2019 +decrease_itm_2020 +decrease_itm_2021 +decrease_itm_2022 +"

R_test_of_significance <- felm(data = dfreg25,weights = dfreg25$population_votingage,formula = as.formula(paste(dep,main_test,Xs_elFE,FEs,cluster,sep = "")))

fstat <- ((sum(R_test_of_significance$residuals^2) - sum(R_pt$residuals^2))/9)/
  (sum(R_pt$residuals^2)/(length(R_pt$residuals)-length(R_pt$coefficients)))

pt_value <- pf(fstat,df2 = 4,length(R_pt$residuals)-length(R_pt$coefficients),lower.tail = FALSE)

# Posttrend test
main_test <- " decrease_itm_2014+decrease_itm_2015 +decrease_itm_2016 + decrease_itm_2017 +"

R_test_of_significance <- felm(data = dfreg25,weights = dfreg25$population_votingage,formula = as.formula(paste(dep,main_test,Xs_elFE,FEs,cluster,sep = "")))

fstat <- ((sum(R_test_of_significance$residuals^2) - sum(R_pt$residuals^2))/9)/
  (sum(R_pt$residuals^2)/(length(R_pt$residuals)-length(R_pt$coefficients)))

pt_value_post <- pf(fstat,df2 = 4,length(R_pt$residuals)-length(R_pt$coefficients),lower.tail = FALSE)


q<-ggplot(data=graphs, aes(x=date, y=coef)) + 
  geom_point(colour = "#B31B1B",size =3) + 
  geom_errorbar(aes(ymin=lb, ymax=up), colour = "#B31B1B", width=100,
                position=position_dodge(0.05))+ 
  labs(x="", y = expression(atop("Coefficients on the interaction between", "year FE and Change.Ded")),
       caption = paste("pre-trend p−value = ",round(pt_value,3),"     -     post-trend p−value = ",round(pt_value_post,3),sep = " ")) + 
  theme_minimal() +
  geom_hline(yintercept=0, linetype=5, color = "#91959C") + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  theme(legend.position = "none") 


# 
if(with_LWP == TRUE){
  dimension <- 2 # highest order term 
  nb_point <- 7# nb on points on CI to search
  source("Rscripts/9bis.LWP.R")
  ggsave(plot = p, file="Output/fig 2.Parallel_trend_lwp.png" , width=8, height=5,bg = "white")
  
}else{
  ggsave(plot = q, file="Output/fig 2.Parallel_trend.png" , width=8, height=5,bg = "white")
  
}

rm(list = setdiff(ls(),c("df","city")))

